/* 2010: Copyright by Massimo Bernaschi, Giorgio Parisi, Leonardo Parisi */
#include "cudamacro.h"
#include "mpiexch.h"

#include "hsgMGPUGM.h"
#include "timing.h"

#include "dictionary.h"
#include "iniparser.h"
#include "util.h"
#if defined(DEBUG)
#include "util/cuPrintf.cu"
#endif

#include "MersenneTwister.h"
#if !defined(HOSTREDU)
#include "threadFenceReduction_kernel.cu"
#endif

extern void StartMpi(int *, int *, int *, char ***);
extern struct mpiexch *SetUpMpi(int, int, int);
extern void ExchMpi(struct mpiexch*);
extern REALENE SumEnergy(REALENE, int);
extern void StopMpi(void);
//-------------------------------------------------------------------------------------------------------
// Usage
//-------------------------------------------------------------------------------------------------------
void Usage(char *program) {
  fprintf(stderr,"Usage: %s -i inputfile [-r spin_file] [-o (over-relax only)] [-b (heat bath only)]\n",program);
}


//-------------------------------------------------------------------------------------------------------
// Init
//-------------------------------------------------------------------------------------------------------
void Init(void){

  if(sizeof(REAL)==sizeof(float)) {
    zero=0.0f;
    one=1.0f;
    two=2.0f;
    three=3.0f;
    zerofive=0.5f;
    zerotwentyfive=0.25f;
  } else if(sizeof(REAL)==sizeof(double)) {
    zero=0.0;
    one=1.0;
    two=2.0;
    three=3.0;
    zerofive=0.5;
    zerotwentyfive=0.25;
  } else {
    writelog(TRUE,APPLICATION_RC,"Invalid size of REAL: %d",sizeof(REAL));
  }

  zeroEne = 0.0; //MESSA SOLO PER IL CALCOLO DELL'ENERGIA IN DOUBLE COME PROVA

  h_Ra = (REAL*) Malloc(ms_s);
  h_Rb = (REAL*) Malloc(ms_s);
  h_Rc = (REAL*) Malloc(ms_s);
  h_Ba = (REAL*) Malloc(ms_s);
  h_Bb = (REAL*) Malloc(ms_s);
  h_Bc = (REAL*) Malloc(ms_s);

  h_Jpx = (REAL*) Malloc(ms_j);
  h_Jpy = (REAL*) Malloc(ms_j);
  h_Jpz = (REAL*) Malloc(ms_j);
  h_Jmx = (REAL*) Malloc(ms_j);
  h_Jmy = (REAL*) Malloc(ms_j);
  h_Jmz = (REAL*) Malloc(ms_j);

  h_Energy = (REAL*) Malloc(ms_s);

  h_sdbL = findlog2(h_side);
  h_sdbLs = findlog2(h_side*h_side);

  h_sdbLm1 = h_sdbL-1;  //log_2(L) -1 used to divide by L/2
  h_sdbLsm1 = h_sdbLs-1;  //log_2(L*L) -1 used to divide by (L*L)/2
  h_maskL = (h_side>>1)-1;
  h_maskL2 = ((h_side*h_side)>>1)-1;
  h_maskL3 = (V(h_side)>>1)-1;
  h_Ldb2 = h_side>>1;
  h_Lsdb2 = (h_side*h_side)>>1;

  h_nspin=h_side*h_side*z_side/2;

  srand48(seed);

}


//-------------------------------------------------------------------------------------------------------
// CopyConstOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyConstOnDevice(void){

  int tmp[2];

  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_side",&h_side,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_nspin",&h_nspin,sizeof(int),0,cudaMemcpyHostToDevice) );

  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbL",&h_sdbL,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLs",&h_sdbLs,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLm1",&h_sdbLm1,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_sdbLsm1",&h_sdbLsm1,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL",&h_maskL,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL2",&h_maskL2,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_maskL3",&h_maskL3,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Ldb2",&h_Ldb2,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Lsdb2",&h_Lsdb2,sizeof(int),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_h",&h_h,sizeof(REAL),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_Beta",&h_Beta,sizeof(REAL),0,cudaMemcpyHostToDevice) );

  tmp[0]=0;
  tmp[1]=-1;

  MY_CUDA_CHECK( cudaMemcpyToSymbol("jindred",tmp,sizeof(tmp),0,cudaMemcpyHostToDevice) );

  tmp[0]=-1;
  tmp[1]=0;

  MY_CUDA_CHECK( cudaMemcpyToSymbol("jindblue",tmp,sizeof(tmp),0,cudaMemcpyHostToDevice) );

  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_one",&one,sizeof(REAL),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_two",&two,sizeof(REAL),0,cudaMemcpyHostToDevice) );
  MY_CUDA_CHECK( cudaMemcpyToSymbol("d_const_zerotwentyfive",&zerotwentyfive,sizeof(REAL),0,cudaMemcpyHostToDevice) );



}


//-------------------------------------------------------------------------------------------------------
// findlog2
//-------------------------------------------------------------------------------------------------------
int findlog2(int n){

  int i=0;

  if(n<1 || n&1) { writelog(TRUE,FINDLOG_RC,"Invalid value %d. Must be > 0 and even\n",n); }

  while(n>0) {
    n>>=1;
    if(n>0) i++;
  }

  return i;
}


//-------------------------------------------------------------------------------------------------------
// nextPow2
//-------------------------------------------------------------------------------------------------------
unsigned int nextPow2( unsigned int x ) {
  --x;
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  return ++x;
}


//-------------------------------------------------------------------------------------------------------
// FreeDevice
//-------------------------------------------------------------------------------------------------------
void FreeDevice(void){

  cudaFree(d_Ra);
  cudaFree(d_Rb);
  cudaFree(d_Rc);
  cudaFree(d_Ba);
  cudaFree(d_Bb);
  cudaFree(d_Bc);

  cudaFree(d_Jpx);
  cudaFree(d_Jpy);
  cudaFree(d_Jpz);
  cudaFree(d_Jmx);
  cudaFree(d_Jmy);
  cudaFree(d_Jmz);

  cudaFree(d_Energy);

}


//-------------------------------------------------------------------------------------------------------
// FreeHost
//-------------------------------------------------------------------------------------------------------
void FreeHost(void){

  Free(h_Ra);
  Free(h_Rb);
  Free(h_Rc);
  Free(h_Ba);
  Free(h_Bb);
  Free(h_Bc);

  Free(h_Jpx);
  Free(h_Jpy);
  Free(h_Jpz);
  Free(h_Jmx);
  Free(h_Jmy);
  Free(h_Jmz);

  Free(h_Energy);

}


//-------------------------------------------------------------------------------------------------------
// Dump_Spin
//-------------------------------------------------------------------------------------------------------
void Dump_Spin(int count){

  FILE* tempfile;
  char namefile[MAXSTRLEN];

  memset(namefile,'\0',MAXSTRLEN); // NON CONTROLLO SE TORNA BENE O MENO

  if(count>=0){

    snprintf(namefile,sizeof(namefile),"spin_dump_%d_%d",count,mpiid);

    MY_CUDA_CHECK( cudaMemcpy( h_Ra, d_Ra+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
    MY_CUDA_CHECK( cudaMemcpy( h_Rb, d_Rb+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
    MY_CUDA_CHECK( cudaMemcpy( h_Rc, d_Rc+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );

    MY_CUDA_CHECK( cudaMemcpy( h_Ba, d_Ba+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
    MY_CUDA_CHECK( cudaMemcpy( h_Bb, d_Bb+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );
    MY_CUDA_CHECK( cudaMemcpy( h_Bc, d_Bc+h_Lsdb2, ms_s, cudaMemcpyDeviceToHost ) );

  } else {
    snprintf(namefile,sizeof(namefile),"spin_dump_start_%d",mpiid);
  }

  tempfile = Fopen(namefile,"w");

  for(int i = 0; i<ns_j; i++){
    if(sizeof(REAL)==sizeof(float)) {
      fprintf(tempfile,"%f %f %f\n%f %f %f\n",
              h_Ra[i],h_Rb[i],h_Rc[i],h_Ba[i],h_Bb[i],h_Bc[i]);
    } else {
      fprintf(tempfile,"%lf %lf %lf\n%lf %lf %lf\n",
              h_Ra[i],h_Rb[i],h_Rc[i],h_Ba[i],h_Bb[i],h_Bc[i]);
    }
  }

  fflush(tempfile);

  fclose(tempfile);

}
//-------------------------------------------------------------------------------------------------------
// Dump_J
//-------------------------------------------------------------------------------------------------------
void Dump_J(void){

  FILE* tempfile;

  tempfile = Fopen("j_dump_start","w");

  if(sizeof(REAL)==sizeof(float)) {

    for(int i = 0; i<ns_j; i++) {
      fprintf(tempfile,"%f %f %f\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);
      fprintf(tempfile,"%f %f %f\n",h_Jmx[i],h_Jmy[i],h_Jmz[i]);
    }

  }else{

    for(int i = 0; i<ns_j; i++) {
      fprintf(tempfile,"%lf %lf %lf\n",h_Jpx[i],h_Jpy[i],h_Jpz[i]);
      fprintf(tempfile,"%lf %lf %lf\n",h_Jmx[i],h_Jmy[i],h_Jmz[i]);
    }

  }


  fflush(tempfile);

  fclose(tempfile);

}


//-------------------------------------------------------------------------------------------------------
// Initialize_Spin
//-------------------------------------------------------------------------------------------------------
void Initialize_Spin(REAL* sa_a, REAL* sa_b, REAL* sa_c){

  REAL a, b, c;
  REAL invnorma, norma;

  for(int i = 0; i<h_nspin; i++){

    norma=two;

    while (norma>one){

      a=(RANDOM01)*two-one;
      b=(RANDOM01)*two-one;
      c=(RANDOM01)*two-one;

      norma=a*a+b*b+c*c;

    }

    invnorma=one/SQRT(norma);

    sa_a[i]=a*invnorma;
    sa_b[i]=b*invnorma;
    sa_c[i]=c*invnorma;

  }

}


//-------------------------------------------------------------------------------------------------------
// Read_Spin
//-------------------------------------------------------------------------------------------------------
void Read_Spin(char *inputfile,
               REAL* sr_a, REAL* sr_b, REAL* sr_c,
               REAL* sb_a, REAL* sb_b, REAL* sb_c) {

  REAL a, b, c;
  int j=0, k=0;
  FILE *fpi;
  fpi=Fopen(inputfile,"r");
  fprintf(stderr,"Reading spin from file %s\n",inputfile);

  for(int i = 0; i< V(h_side); i++){
    if(sizeof(REAL)==sizeof(float)) {
      fscanf(fpi,"%f %f %f",&a,&b,&c);
    } else {
      fscanf(fpi,"%lf %lf %lf",&a,&b,&c);
    }
    if(i&1) {
      sb_a[j]=a;
      sb_b[j]=b;
      sb_c[j]=c;
      j++;
    } else {
      sr_a[k]=a;
      sr_b[k]=b;
      sr_c[k]=c;
      k++;
    }
  }
  fclose(fpi);

}


//-------------------------------------------------------------------------------------------------------
// CopyDataOnDevice
//-------------------------------------------------------------------------------------------------------
void CopyDataOnDevice(void){

  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Ra, ms_s_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Rb, ms_s_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Rc, ms_s_mpi ) );

  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Ba, ms_s_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Bb, ms_s_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Bc, ms_s_mpi ) );

  MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_Lsdb2, h_Ra, ms_s, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_Lsdb2, h_Rb, ms_s, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_Lsdb2, h_Rc, ms_s, cudaMemcpyHostToDevice ) );

  MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_Lsdb2, h_Ba, ms_s, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_Lsdb2, h_Bb, ms_s, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_Lsdb2, h_Bc, ms_s, cudaMemcpyHostToDevice ) );

  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpx, ms_j_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpy, ms_j_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jpz, ms_j_mpi ) );

  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmx, ms_j_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmy, ms_j_mpi ) );
  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Jmz, ms_j_mpi ) );

  MY_CUDA_CHECK( cudaMemcpy( d_Jpx+h_Lsdb2, h_Jpx, ms_j, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jpy+h_Lsdb2, h_Jpy, ms_j, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jpz+h_Lsdb2, h_Jpz, ms_j, cudaMemcpyHostToDevice ) );

  MY_CUDA_CHECK( cudaMemcpy( d_Jmx+h_Lsdb2, h_Jmx, ms_j, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmy+h_Lsdb2, h_Jmy, ms_j, cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmz+h_Lsdb2, h_Jmz, ms_j, cudaMemcpyHostToDevice ) );

  MY_CUDA_CHECK( cudaMalloc( (void **) &d_Energy, ms_s ) );
  MY_CUDA_CHECK( cudaMemset( d_Energy, 0, ms_s ) );
}


//-------------------------------------------------------------------------------------------------------
// Initialize_J
//-------------------------------------------------------------------------------------------------------
void Initialize_J(void){

  int j;

  switch (init_type_j){

  case 0:
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jpy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jpz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jmx[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jmy[j] = ((RANDOM01)>zerofive) ? -1 : 1;
      h_Jmz[j] = ((RANDOM01)>zerofive) ? -1 : 1;
    }

    break;
  case 1:
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = h_Jpy[j] = h_Jpz[j] =
        h_Jmx[j] = h_Jmy[j] = h_Jmz[j] = one;
    }
    break;
  case -1:
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = h_Jpy[j] = h_Jpz[j] =
        h_Jmx[j] = h_Jmy[j] = h_Jmz[j] = -one;
    }

    break;
  case 2: //Caso Gaussiano
    for (j=0;j<ns_j;j++) {
      h_Jpx[j] = rand_gauss();
      h_Jpy[j] = rand_gauss();
      h_Jpz[j] = rand_gauss();
      h_Jmx[j] = rand_gauss();
      h_Jmy[j] = rand_gauss();
      h_Jmz[j] = rand_gauss();
    }
    break;
  default: {
    writelog(TRUE,APPLICATION_RC,"Invalid initialization for J: %d", init_type_j);
    break;
  }
  }

}


//-------------------------------------------------------------------------------------------------------
// rand_gauss
//-------------------------------------------------------------------------------------------------------
REAL rand_gauss(void){

  REAL r,cosr;
  REAL a=one-RANDOM01;
  if(a<EPSILON) { a=EPSILON; }

  r=SQRT(-2.*LOG(a));

  cosr=COS((2*M_PI)*(RANDOM01));

  return r*cosr;

}


//-------------------------------------------------------------------------------------------------------
// Parameters_Reading
//-------------------------------------------------------------------------------------------------------
void Parameters_Reading(char *inputfile){

  dictionary *ini;
  char key[MAXSTRLEN];

  if(inputfile==NULL) { writelog(TRUE,APPLICATION_RC,"No inputfile"); }

  ini = iniparser_load(inputfile);

  if(ini==NULL) { writelog(TRUE,APPLICATION_RC,"Cannot parse file: %s\n", inputfile); }

  READINTFI(h_side,"Side");
  READINTFI(init_type_j,"Initialization type for J");
  READINTFI(nsteps,"Number of steps");
  READINTFI(overrelaxiter,"Number of overrelaxation iterations");
  READINTFI(outfreq,"Output frequency");
  READINTFI(dumpfreq,"Dump configuration frequency");
  READREALFI(h_h,"Magnetic field");
  READREALFI(h_Beta,"Beta (1/T)");
  READINTFI(seed,"Random seed");

  DEVICEREADINTFI(numBlocks,"Number of blocks");
  DEVICEREADINTFI(numThreadsPerBlock,"Number of threads for block");
  DEVICEREADINTFI(whichgpu,"Which GPU");

  DEVICEREADINTFI(RandNumBlocks,"RandGene num of blocks");
  DEVICEREADINTFI(RandNumThreadsPerBlock,"RandGene num of threads for block");

}


//-------------------------------------------------------------------------------------------------------
// Calculate_Sizes
//-------------------------------------------------------------------------------------------------------
void Calculate_Sizes(void){

  ns_j = h_side*h_side*z_side/2;

  ms_s = (h_side*h_side*z_side/2)*sizeof(REAL);
  ms_s_mpi = ((h_side*h_side*z_side/2)+(h_side*h_side))*sizeof(REAL);

  ms_j = ns_j*sizeof(REAL);
  ms_j_mpi = ((h_side*h_side*z_side/2)+(h_side*h_side))*sizeof(REAL);

}


//-------------------------------------------------------------------------------------------------------
// getNumBlocksAndThreads
//-------------------------------------------------------------------------------------------------------
void getNumBlocksAndThreads(int n, int maxBlocks, int maxThreads, int &blocks, int &threads) {

  if (n == 1){
    threads = 1;
    blocks = 1;
  } else{
    threads = (n < maxThreads*2) ? nextPow2(n / 2) : maxThreads;
    blocks = max(1, n / (threads * 2));
  }

  blocks = min(maxBlocks, blocks);
}


//-------------------------------------------------------------------------------------------------------
// main
//-------------------------------------------------------------------------------------------------------
int main(int argc, char **argv){

  char *inputfile = NULL;
  char *po;
  int istep;
  int nblocks_red, nthreads_red;
  int cudadevice, ngpu;
  double ttimeover=0., ttimehb=0.;
  REAL normene;
  int j;
  char *readspin=NULL;
#if !defined(HOSTREDU)
  unsigned int retCnt;
#endif
  int myid=0, nprocs=1;
  int lpoff;
  struct mpiexch *dataexch;

  TIMER_DEF;

  for(int i = 1; i < argc; i++) {
    po = argv[i];
    if (*po++ == '-') {
      switch (*po++) {
      case '?':
        SKIPBLANK
          Usage(argv[0]);
        exit(OK);
        break;
      case 'v':
        SKIPBLANK
          verbose=TRUE;
        break;
      case 'i':
        SKIPBLANK
          inputfile=Strdup(po);
        break;
      case 'o':
        SKIPBLANK
          onlyOverrelaxation=true;
        break;
      case 'r':
        SKIPBLANK
          readspin=Strdup(po);
        break;

      case 'b':
        SKIPBLANK
          onlyHeatBath=true;
        break;
      default:
        Usage(argv[0]);
        exit(OK);
        break;
      }
    }
  }

  StartMpi(&myid, &nprocs,&argc,&argv);
  mpiid=myid;
#if defined(DEBUG)
  cudaPrintfInit();
#endif
  cudaGetDeviceCount(&ngpu);
  if(ngpu<1) {
    writelog(TRUE,APPLICATION_RC,"No gpu available! %d\n",ngpu);
  }
  Parameters_Reading(inputfile);

  seed=(seed+myid*nprocs)%MAXSEED;

  dataexch=SetUpMpi(myid,nprocs,h_side);
  z_side=h_side/nprocs;
  if((z_side*nprocs)<h_side) {
    writelog(TRUE,APPLICATION_RC,"Invalid number of processors (%d) or system size (%d). Number of processors must be an exact divider of the system size",
                                  nprocs, h_side);
  }

  if(whichgpu>(ngpu-1)) {
        writelog(TRUE,APPLICATION_RC,
        "invalid number of gpu: %d must be > 0 and <= %d\n",ngpu);
  }

  if(whichgpu>0) {
        printf("There are %d GPUs!, forcing usage of n. <%d>\n",ngpu,whichgpu);
        cudaSetDevice(whichgpu);
  }

  if (h_side<=0 || !pow_of_two(h_side)) {
    writelog(TRUE,APPLICATION_RC,"Invalid side: must be > 0 and power of two\n");
  }

  Calculate_Sizes();
  Init();

  lpoff=(z_side-1)*h_Lsdb2;

  if(readspin==NULL) {

    Initialize_Spin(h_Ra,h_Rb,h_Rc);
    Initialize_Spin(h_Ba,h_Bb,h_Bc);
    Dump_Spin(-1);

  } else {
    Read_Spin(readspin,h_Ra,h_Rb,h_Rc, h_Ba,h_Bb,h_Bc);
  }

  Initialize_J();

  // Dump_J();

  initRand(V(h_side));

  CopyConstOnDevice();

  CopyDataOnDevice();

  getNumBlocksAndThreads(h_nspin, numBlocks, numThreadsPerBlock,nblocks_red, nthreads_red);

  normene=one/(three*V(h_side));

  h_Energy[0]=zero;

  ComputeHEnergy(h_Ra, h_Rb, h_Rc, h_Ba, h_Bb, h_Bc,h_Jpx, h_Jpy, h_Jpz, h_Jmx, h_Jmy, h_Jmz, h_Energy,TRUE);
  h_Energy[0]=SumEnergy(h_Energy[0], myid);
  if(myid==0 && nprocs==1) {
    printf("Initial Energy per spin on host: %f\n", normene*h_Energy[0]);
  }

  memcpy(dataexch->downsend,h_Ra,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+h_Lsdb2,h_Rb,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(2*h_Lsdb2),h_Rc,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(3*h_Lsdb2),h_Ba,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(4*h_Lsdb2),h_Bb,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(5*h_Lsdb2),h_Bc,h_Lsdb2*sizeof(REAL));

  memcpy(dataexch->upsend,h_Ra+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+h_Lsdb2,h_Rb+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(2*h_Lsdb2),h_Rc+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(3*h_Lsdb2),h_Ba+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(4*h_Lsdb2),h_Bb+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(5*h_Lsdb2),h_Bc+lpoff,h_Lsdb2*sizeof(REAL));


  dataexch->nbyte=2*h_Lsdb2*sizeof(REAL)*NCOMP;
  dataexch->post=FALSE;
  ExchMpi(dataexch);

  MY_CUDA_CHECK( cudaMemcpy( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Rb, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Rc, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Ba, dataexch->downrecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Bb, dataexch->downrecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Bc, dataexch->downrecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );

  MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv,
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2,
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv+(3*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+(4*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_nspin+h_Lsdb2, dataexch->uprecv+(5*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );

  memcpy(dataexch->downsend,h_Jpx,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+h_Lsdb2,h_Jpy,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(2*h_Lsdb2),h_Jpz,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(3*h_Lsdb2),h_Jmx,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(4*h_Lsdb2),h_Jmy,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->downsend+(5*h_Lsdb2),h_Jmz,h_Lsdb2*sizeof(REAL));

  memcpy(dataexch->upsend,h_Jpx+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+h_Lsdb2,h_Jpy+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(2*h_Lsdb2),h_Jpz+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(3*h_Lsdb2),h_Jmx+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(4*h_Lsdb2),h_Jmy+lpoff,h_Lsdb2*sizeof(REAL));
  memcpy(dataexch->upsend+(5*h_Lsdb2),h_Jmz+lpoff,h_Lsdb2*sizeof(REAL));

  dataexch->post=FALSE;
  ExchMpi(dataexch);

  MY_CUDA_CHECK( cudaMemcpy( d_Jpx, dataexch->downrecv, h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jpy, dataexch->downrecv+h_Lsdb2, h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jpz, dataexch->downrecv+(2*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmx, dataexch->downrecv+(3*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmy, dataexch->downrecv+(4*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmz, dataexch->downrecv+(5*h_Lsdb2), h_Lsdb2*sizeof(REAL),
                             cudaMemcpyHostToDevice ) );

  MY_CUDA_CHECK( cudaMemcpy( d_Jpx+h_nspin+h_Lsdb2, dataexch->uprecv,
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jpy+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2,
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jpz+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmx+h_nspin+h_Lsdb2, dataexch->uprecv+(3*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmy+h_nspin+h_Lsdb2, dataexch->uprecv+(4*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
  MY_CUDA_CHECK( cudaMemcpy( d_Jmz+h_nspin+h_Lsdb2, dataexch->uprecv+(5*h_Lsdb2),
                             h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );

  dataexch->post=FALSE;
  dataexch->nbyte=h_Lsdb2*sizeof(REAL)*NCOMP;

  MY_CUDA_CHECK( cudaMalloc((void **)&r_odata,numBlocks*sizeof(REAL)) );

  dim3 dimGrid(numBlocks,1,1);
  dim3 dimBlock(numThreadsPerBlock,1,1);

  ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                        d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                        d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                        d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                        d_Energy, TRUE);
#if defined(HOSTREDU)
  MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );

  h_totene=zeroEne;

  for(int i=0; i<h_nspin; i++) { h_totene+=h_Energy[i]; }

#else
  MY_CUDA_CHECK( cudaThreadSynchronize() );
  retCnt=0;
  cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
  reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
  MY_CUDA_CHECK( cudaThreadSynchronize() );
  MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );
#endif
  h_totene=SumEnergy(h_totene, myid);
  if(myid==0) {
    printf("Initial Energy per spin: %f\n", normene*h_totene);
  }

#if defined(FERMI)
  cudaFuncSetCacheConfig(OverrelaxRed, cudaFuncCachePreferL1 );
  cudaFuncSetCacheConfig(OverrelaxBlue, cudaFuncCachePreferL1 );
#endif

  for(istep=0; istep<nsteps; istep++) {

    if(!onlyHeatBath){

      TIMER_START;

      for(j=0; j<overrelaxiter; j++) {

#if defined(SINGLEOVERKERNEL)
        Overrelaxation<<< dimGrid,dimBlock>>>(d_Ra, d_Rb, d_Rc,
                                              d_Ba, d_Bb, d_Bc,
                                              d_Jpx, d_Jpy, d_Jpz,
                                              d_Jmx, d_Jmy, d_Jmz, TRUE);

        MY_CUDA_CHECK( cudaThreadSynchronize() );

        Overrelaxation<<< dimGrid,dimBlock>>>(d_Ba, d_Bb, d_Bc,
                                              d_Ra, d_Rb, d_Rc,
                                              d_Jmx, d_Jmy, d_Jmz,
                                              d_Jpx, d_Jpy, d_Jpz, FALSE);

        MY_CUDA_CHECK( cudaThreadSynchronize() );
#else
        OverrelaxRed<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                             d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                             d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                             d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2);
        //        MY_CUDA_CHECK( cudaThreadSynchronize() );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->downsend, d_Ra+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->downsend+h_Lsdb2, d_Rb+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->downsend+(2*h_Lsdb2), d_Rc+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );

        MY_CUDA_CHECK( cudaMemcpy( dataexch->upsend, d_Ra+h_nspin, h_Lsdb2*sizeof(REAL),
                                   cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->upsend+h_Lsdb2, d_Rb+h_nspin,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->upsend+(2*h_Lsdb2), d_Rc+h_nspin,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );

        ExchMpi(dataexch);

        MY_CUDA_CHECK( cudaMemcpy( d_Ra, dataexch->downrecv, h_Lsdb2*sizeof(REAL),
                                   cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Rb, dataexch->downrecv+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Rc, dataexch->downrecv+(2*h_Lsdb2),
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );

        MY_CUDA_CHECK( cudaMemcpy( d_Ra+h_nspin+h_Lsdb2, dataexch->uprecv,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Rb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Rc+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2),
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );

        OverrelaxBlue<<< dimGrid,dimBlock>>>(d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                             d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                             d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                             d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2);

        MY_CUDA_CHECK( cudaThreadSynchronize() );

        MY_CUDA_CHECK( cudaMemcpy( dataexch->downsend, d_Ba+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->downsend+h_Lsdb2, d_Bb+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->downsend+(2*h_Lsdb2), d_Bc+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );

        MY_CUDA_CHECK( cudaMemcpy( dataexch->upsend, d_Ba+h_nspin,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->upsend+h_Lsdb2, d_Bb+h_nspin,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );
        MY_CUDA_CHECK( cudaMemcpy( dataexch->upsend+(2*h_Lsdb2), d_Bc+h_nspin,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyDeviceToHost ) );

        ExchMpi(dataexch);

        MY_CUDA_CHECK( cudaMemcpy( d_Ba, dataexch->downrecv,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Bb, dataexch->downrecv+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Bc, dataexch->downrecv+(2*h_Lsdb2),
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );

        MY_CUDA_CHECK( cudaMemcpy( d_Ba+h_nspin+h_Lsdb2, dataexch->uprecv,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Bb+h_nspin+h_Lsdb2, dataexch->uprecv+h_Lsdb2,
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
        MY_CUDA_CHECK( cudaMemcpy( d_Bc+h_nspin+h_Lsdb2, dataexch->uprecv+(2*h_Lsdb2),
                                   h_Lsdb2*sizeof(REAL), cudaMemcpyHostToDevice ) );
#endif
      }

      TIMER_STOP;
      ttimeover+=TIMER_ELAPSED;
    }

    if(!onlyOverrelaxation) {

      TIMER_START;

      /* Here it goes the heat-bath */
      /* first of all, random number generation */
      /* then the heat-bath itself */
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ONEONE,d_cRand);
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_aRand);
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROFIVE,d_bRand);
      ComputeRandomGPU(RandNumBlocks,RandNumThreadsPerBlock,ZEROONE,d_Rand);

      if(istep==0&&FALSE){
        printRand(randNum,d_aRand,"daRand3");
        printRand(randNum,d_bRand,"dbRand3");
        printRand(randNum,d_cRand,"dcRand3");
        printRand(randNum,d_Rand,"dRand3");
        exit(1);
      }


#if defined(SINGLEHBKERNEL)
      HeatBath<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,
                                       d_Ba, d_Bb, d_Bc,
                                       d_Jpx, d_Jpy, d_Jpz,
                                       d_Jmx, d_Jmy, d_Jmz,
                                       d_aRand,d_bRand,d_cRand,
                                       d_Rand,TRUE);

      MY_CUDA_CHECK( cudaThreadSynchronize() );

      HeatBath<<< dimGrid,dimBlock >>>(d_Ba, d_Bb, d_Bc,
                                       d_Ra, d_Rb, d_Rc,
                                       d_Jmx, d_Jmy, d_Jmz,
                                       d_Jpx, d_Jpy, d_Jpz,
                                       d_aRand+h_nspin,d_bRand+h_nspin,
                                       d_cRand+h_nspin,
                                       d_Rand+h_nspin,FALSE);

      MY_CUDA_CHECK( cudaThreadSynchronize() );
#else
      HeatBathRed<<< dimGrid,dimBlock >>>(d_Ra, d_Rb, d_Rc,
                                          d_Ba, d_Bb, d_Bc,
                                          d_Jpx, d_Jpy, d_Jpz,
                                          d_Jmx, d_Jmy, d_Jmz,
                                          d_aRand,d_bRand,d_cRand,
                                          d_Rand);

      MY_CUDA_CHECK( cudaThreadSynchronize() );

      HeatBathBlue<<< dimGrid,dimBlock >>>(d_Ba, d_Bb, d_Bc,
                                           d_Ra, d_Rb, d_Rc,
                                           d_Jmx, d_Jmy, d_Jmz,
                                           d_Jpx, d_Jpy, d_Jpz,
                                           d_aRand+h_nspin,d_bRand+h_nspin,
                                           d_cRand+h_nspin,
                                           d_Rand+h_nspin);

      MY_CUDA_CHECK( cudaThreadSynchronize() );

#endif

      TIMER_STOP;
      ttimehb+=TIMER_ELAPSED;

    }

    if(istep && (istep%outfreq)==0) {
      ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                        d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                        d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                        d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                        d_Energy, TRUE);

      MY_CUDA_CHECK( cudaThreadSynchronize() );


#if defined(HOSTREDU)

      MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy,
                                 h_nspin*sizeof(REAL),
                                 cudaMemcpyDeviceToHost) );
      h_totene=zeroEne;

      for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];

#else
      retCnt=0;
      cudaMemcpyToSymbol("retirementCount", &retCnt, sizeof(unsigned int), 0, cudaMemcpyHostToDevice);
      reduceSinglePass(h_nspin, nthreads_red, nblocks_red, d_Energy, r_odata);
      MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );

#endif
      h_totene=SumEnergy(h_totene, myid);
      if(myid==0) {
        printf("Energy per spin at step %d: %f\n",istep,normene*h_totene);
      }
      //place holder for energy evaluation
    }

    if(istep && (istep%dumpfreq)==0) {
      //Dump_Spin(istep);
    }

  }

  if(myid==0) {
    if(!onlyHeatBath){
      printf("Total time for %d iterations of overrelaxation on a %d^3 cube: %f seconds\n", nsteps*overrelaxiter, h_side, ttimeover*0.000001);

      printf("Time for spin: %f microseconds\n",ttimeover/(overrelaxiter*(double)nsteps*V(h_side)));
    }

    if(!onlyOverrelaxation){
      printf("Total time for %d iterations of heat bath on a %d^3 cube: %f seconds\n", nsteps, h_side, ttimehb*0.000001);

      printf("Time for spin: %f microseconds\n",ttimehb/((double)nsteps*V(h_side)));
    }
  }

  ComputeEnergy<<< dimGrid,dimBlock >>>(d_Ra+h_Lsdb2, d_Rb+h_Lsdb2, d_Rc+h_Lsdb2,
                                        d_Ba+h_Lsdb2, d_Bb+h_Lsdb2, d_Bc+h_Lsdb2,
                                        d_Jpx+h_Lsdb2, d_Jpy+h_Lsdb2, d_Jpz+h_Lsdb2,
                                        d_Jmx+h_Lsdb2, d_Jmy+h_Lsdb2, d_Jmz+h_Lsdb2,
                                        d_Energy, TRUE);

#if defined(HOSTREDU)

  MY_CUDA_CHECK( cudaMemcpy(h_Energy, d_Energy, h_nspin*sizeof(REAL),cudaMemcpyDeviceToHost) );

  h_totene=zero;

  for(int i=0; i<h_nspin; i++) h_totene+=h_Energy[i];

#else

  MY_CUDA_CHECK( cudaThreadSynchronize() );
  reduceSinglePass(ns_j, nthreads_red, nblocks_red,d_Energy, r_odata);
  MY_CUDA_CHECK( cudaMemcpy(&h_totene, r_odata, sizeof(REAL),cudaMemcpyDeviceToHost) );

#endif
  h_totene=SumEnergy(h_totene, myid);
  if(myid==0) {
    printf("Final Energy per spin: %f\n",normene*h_totene);
  }

  Dump_Spin(nsteps);

  FreeDevice();

  FreeHost();

  FreeRand();
#if defined(DEBUG)
  cudaPrintfDisplay();
  cudaPrintfEnd();
#endif
  StopMpi();

  return 0;

}


//-------------------------------------------------------------------------------------------------------
// Overrelaxation
//-------------------------------------------------------------------------------------------------------
#if defined(SINGLEOVERKERNEL)
__global__ void Overrelaxation(REAL *a, REAL *b, REAL *c,
                               REAL *an, REAL *bn, REAL *cn,
                               REAL *Jpx, REAL *Jpy, REAL *Jpz,
                               REAL *Jmx, REAL *Jmy, REAL *Jmz,
                               int even) {

  unsigned int iy, iz;
  int diff;
  REAL as, bs, cs;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x+d_Lsdb2; /*start on second plane (by using d_Lsdb2)
                                          since first plane is for the boundary */


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(((iz&1)^(iy&1))^even);

    spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) & d_maskL3;
    smz=(i-d_Lsdb2) & d_maskL3;

    as=an[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
      an[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
      an[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
      an[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
      an[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
      an[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];

    bs=bn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
      bn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
      bn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
      bn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
      bn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
      bn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];

    cs=cn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
      cn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
      cn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
      cn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
      cn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
      cn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];

    cs+=d_h;

    //Overrelaxation

    factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);

    a[i]=as*factor-a[i];
    b[i]=bs*factor-b[i];
    c[i]=cs*factor-c[i];

  }

}

#else

//-------------------------------------------------------------------------------------------------------
// OverrelaxRed
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxRed(REAL *a, REAL *b, REAL *c,
                             REAL *an, REAL *bn, REAL *cn,
                             REAL *Jpx, REAL *Jpy, REAL *Jpz,
                             REAL *Jmx, REAL *Jmy, REAL *Jmz) {

  unsigned int iy, iz;
  int diff;
  REAL as, bs, cs;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
                                          since first plane is for the boundary */


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(iz&1)^(iy&1);

    spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) /* & d_maskL3 */;
    smz=(i-d_Lsdb2) /* & d_maskL3 */;

    as=an[spx]*Jpx[i]+
      an[spy]*Jpy[i]+
      an[spz]*Jpz[i]+
      an[smx]*Jmx[i]+
      an[smy]*Jmy[i]+
      an[smz]*Jmz[i];

    bs=bn[spx]*Jpx[i]+
      bn[spy]*Jpy[i]+
      bn[spz]*Jpz[i]+
      bn[smx]*Jmx[i]+
      bn[smy]*Jmy[i]+
      bn[smz]*Jmz[i];

    cs=cn[spx]*Jpx[i]+
      cn[spy]*Jpy[i]+
      cn[spz]*Jpz[i]+
      cn[smx]*Jmx[i]+
      cn[smy]*Jmy[i]+
      cn[smz]*Jmz[i];

    cs+=d_h;

    //Overrelaxation

    factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/(as*as+bs*bs+cs*cs);

#if defined(DEBUG)
    printf("Reds [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
#endif

    a[i]=as*factor-a[i];
    b[i]=bs*factor-b[i];
    c[i]=cs*factor-c[i];

#if defined(DEBUG)
    printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif


  }

}


//-------------------------------------------------------------------------------------------------------
// OverrelaxBlue
//-------------------------------------------------------------------------------------------------------
__global__ void OverrelaxBlue(REAL *a, REAL *b, REAL *c,
                              REAL *an, REAL *bn, REAL *cn,
                              REAL *Jpx, REAL *Jpy, REAL *Jpz,
                              REAL *Jmx, REAL *Jmy, REAL *Jmz) {

  unsigned int iy, iz;
  int diff;
  REAL as, bs, cs;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
                                          since first plane is for the boundary */


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(iz&1)^(iy&1);

    spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) /* & d_maskL3 */;
    smz=(i-d_Lsdb2) /* & d_maskL3 */;

    as=an[spx]*Jpx[spx]+
      an[spy]*Jpy[spy]+
      an[spz]*Jpz[spz]+
      an[smx]*Jmx[smx]+
      an[smy]*Jmy[smy]+
      an[smz]*Jmz[smz];

    bs=bn[spx]*Jpx[spx]+
      bn[spy]*Jpy[spy]+
      bn[spz]*Jpz[spz]+
      bn[smx]*Jmx[smx]+
      bn[smy]*Jmy[smy]+
      bn[smz]*Jmz[smz];

    cs=cn[spx]*Jpx[spx]+
      cn[spy]*Jpy[spy]+
      cn[spz]*Jpz[spz]+
      cn[smx]*Jmx[smx]+
      cn[smy]*Jmy[smy]+
      cn[smz]*Jmz[smz];

    cs+=d_h;

    //Overrelaxation

    factor=d_const_two*(as*a[i]+bs*b[i]+cs*c[i])/
      (as*as+bs*bs+cs*cs);
#if defined(DEBUG)
    printf("Blues [%d]: %7.5f %7.5f %7.5f, xS= %7.5f %7.5f %7.5f, factor: %7.5f ",i,a[i],b[i],c[i], as,bs,cs, factor);
    cuPrintf("i=%d, spx=%d, spy=%d, spz=%d, smx=%d, smy=%d, smz=%d, factor=%f", i, spx, spy, spz, smx, smy, smz, factor);
    cuPrintf("a=%f, b=%f, c=%f, anspx=%f, anspy=%f, anspz=%f, ansmx=%f, ansmy=%f, ansmz=%f, ", a[i], b[i], c[i], an[spx],an[spy],an[spz], an[smx],an[smy],an[smz]);

    cuPrintf("bnspx=%f, bnspy=%f, bnspz=%f, bnsmx=%f, bnsmy=%f, bnsmz=%f, ", bn[spx],bn[spy],bn[spz], bn[smx],bn[smy],bn[smz]);
    cuPrintf("cnspx=%f, cnspy=%f, cnspz=%f, cnsmx=%f, cnsmy=%f, cnsmz=%f, ", cn[spx],cn[spy],cn[spz], cn[smx],cn[smy],cn[smz]);
    cuPrintf("Jpx=%f, Jpy=%f, Jpz=%f, Jmx=%f, Jmy=%f, Jmz=%f\n" , Jpx[spx],Jpy[spy],Jpz[spz], Jmx[smx],Jmy[smy],Jmz[smz]);
#endif
    a[i]=as*factor-a[i];
    b[i]=bs*factor-b[i];
    c[i]=cs*factor-c[i];
#if defined(DEBUG)
    printf("out: %7.5f %7.5f\n",a[i],b[i],c[i]);
#endif

  }

}

#endif


//-------------------------------------------------------------------------------------------------------
// ComputeEnergy
//-------------------------------------------------------------------------------------------------------
__global__ void ComputeEnergy(REAL *a, REAL *b, REAL *c,
                              REAL *an, REAL *bn, REAL *cn,
                              REAL *Jpx, REAL *Jpy, REAL *Jpz,
                              REAL *Jmx, REAL *Jmy, REAL *Jmz,
                              REAL *energy,
                              int even) {

  unsigned int iy, iz;
  int diff;
  int spx, smx, spy, smy, spz, smz;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x; /*start on second plane (by using d_Lsdb2)
                                          since first plane is for the boundary */


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(((iz&1)^(iy&1))^even);

    spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) /* & d_maskL3 */;
    smz=(i-d_Lsdb2) /* & d_maskL3 */;


    energy[i]=a[i]*(an[spx]*Jpx[i]+
                    an[spy]*Jpy[i]+
                    an[spz]*Jpz[i]+
                    an[smx]*Jmx[i]+
                    an[smy]*Jmy[i]+
                    an[smz]*Jmz[i])+
      b[i]*(bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i])+
      c[i]*(cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i]);
#if 0
    cuPrintf("i=%d, spx=%d, spy=%d, spz=%d, smx=%d, smy=%d, smz=%d, ", i, spx, spy, spz, smx, smy, smz);
    cuPrintf("a=%f, b=%f, c=%f, anspx=%f, anspy=%f, anspz=%f, ansmx=%f, ansmy=%f, ansmz=%f, ", a[i], b[i], c[i], an[spx],an[spy],an[spz], an[smx],an[smy],an[smz]);

    cuPrintf("bnspx=%f, bnspy=%f, bnspz=%f, bnsmx=%f, bnsmy=%f, bnsmz=%f, ", bn[spx],bn[spy],bn[spz], bn[smx],bn[smy],bn[smz]);
    cuPrintf("cnspx=%f, cnspy=%f, cnspz=%f, cnsmx=%f, cnsmy=%f, cnsmz=%f, ", cn[spx],cn[spy],cn[spz], cn[smx],cn[smy],cn[smz]);
    cuPrintf("Jpx=%f, Jpy=%f, Jpz=%f, Jmx=%f, Jmy=%f, Jmz=%f\n" , Jpx[i],Jpy[i],Jpz[i], Jmx[i],Jmy[smy],Jmz[i]);
#endif
  }

}


//-------------------------------------------------------------------------------------------------------
// ComputeHEnergy
//-------------------------------------------------------------------------------------------------------
void ComputeHEnergy(REAL *a, REAL *b, REAL *c,
                    REAL *an, REAL *bn, REAL *cn,
                    REAL *Jpx, REAL *Jpy, REAL *Jpz,
                    REAL *Jmx, REAL *Jmy, REAL *Jmz,
                    REAL *energy,
                    int even) {

  unsigned int iy, iz;
  int diff;
  int spx, smx, spy, smy, spz, smz;
  REAL e;

  int i;

  for(i=0; i<h_nspin; i++) {
    iz=i>>h_sdbLsm1;
    iy=(i - (iz<<h_sdbLsm1))>>h_sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(((iz&1)^(iy&1))^even);

    spx=(i >> (diff-1)) + (((i-(i&h_maskL)) + ((i+1)&h_maskL)) >> (-diff));
    smx=(((i-(i&h_maskL)) + ((i-1)&h_maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&h_maskL2))+((i+h_Ldb2)&h_maskL2);
    smy=(i-(i&h_maskL2))+((i-h_Ldb2)&h_maskL2);
    spz=(i+h_Lsdb2) & h_maskL3;
    smz=(i-h_Lsdb2) & h_maskL3;

    e=a[i]*(an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i])+
      b[i]*(bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i])+
      c[i]*(cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i]);
    energy[0]+=e;

  }

}


//-------------------------------------------------------------------------------------------------------
// ComputeRandomGPU
//-------------------------------------------------------------------------------------------------------
void ComputeRandomGPU(int numBlocksRand, int numThreadsPerBlockRand, int mode, REAL* array){

  seedRand(0);

  MY_CUDA_CHECK( cudaThreadSynchronize() );

  dim3 dimGrid(numBlocksRand,1,1);

  dim3 dimBlock(numThreadsPerBlockRand,1,1);

  RandomGPU<<<dimGrid, numThreadsPerBlockRand>>>(array,mode);

  //cutilCheckMsg("RandomGPU() execution failed\n");

  MY_CUDA_CHECK( cudaThreadSynchronize() );

}


//-------------------------------------------------------------------------------------------------------
// FreeRand
//-------------------------------------------------------------------------------------------------------
void FreeRand(void){

  MY_CUDA_CHECK( cudaFree(d_aRand) );
  MY_CUDA_CHECK( cudaFree(d_bRand) );
  MY_CUDA_CHECK( cudaFree(d_cRand) );
  MY_CUDA_CHECK( cudaFree(d_Rand) );

}


//-------------------------------------------------------------------------------------------------------
// initRand
//-------------------------------------------------------------------------------------------------------
int initRand(int randomNumber){


  int n_per_rng;
  // int randNum; DA DISCOMMENTARE DOPO PROVA MESSO GLOBAL

  n_per_rng = iAlignUp(iDivUp(randomNumber, MT_RNG_COUNT), 2);

  randNum = MT_RNG_COUNT * n_per_rng;

  MY_CUDA_CHECK( cudaMalloc((void **)&d_aRand, randNum * sizeof(REAL)) );
  MY_CUDA_CHECK( cudaMalloc((void **)&d_bRand, randNum * sizeof(REAL)) );
  MY_CUDA_CHECK( cudaMalloc((void **)&d_cRand, randNum * sizeof(REAL)) );
  MY_CUDA_CHECK( cudaMalloc((void **)&d_Rand, randNum * sizeof(REAL)) );

  MY_CUDA_CHECK( cudaMemcpyToSymbol("NPerRng",&n_per_rng,sizeof(int),0,cudaMemcpyHostToDevice) );

  return randNum;

}


//-------------------------------------------------------------------------------------------------------
// seedRand
//-------------------------------------------------------------------------------------------------------
void seedRand(int seed){

  if(seed==0) seed=(int)(drand48()*1023); // UNA SCHIFEZZA

  static mt_struct_stripped h_MT[MT_RNG_COUNT];
  static mt_struct MT[MT_RNG_COUNT];

  FILE *fd = Fopen("./data/MersenneTwister.raw", "rb");

  for (int i = 0; i < MT_RNG_COUNT; i++){
    //Inline structure size for compatibility,
    //since pointer types are 8-byte on 64-bit systems (unused *state variable)
    Fread(MT + i, 16 /* sizeof(mt_struct) */ * sizeof(int), 1, fd);

  }

  fclose(fd);

  fd = Fopen("./data/MersenneTwister.dat", "rb");

  Fread(h_MT, sizeof(h_MT), 1, fd);

  fclose(fd);

  //Need to be thread-safe
  mt_struct_stripped *MT_temp = (mt_struct_stripped *)Malloc(MT_RNG_COUNT * sizeof(mt_struct_stripped));

  for(int i = 0; i < MT_RNG_COUNT; i++){
    MT_temp[i] = h_MT[i];
    MT_temp[i].seed = seed;
  }

  MY_CUDA_CHECK( cudaMemcpyToSymbol("ds_MT",MT_temp,sizeof(h_MT),0,cudaMemcpyHostToDevice) );

  free(MT_temp);

}


//-------------------------------------------------------------------------------------------------------
// printRand
//-------------------------------------------------------------------------------------------------------
void printRand(int randomNumber,REAL* array,char* namefile){

  REAL * h_Rand  = (REAL *)malloc(randomNumber * sizeof(REAL));

  MY_CUDA_CHECK( cudaMemcpy(h_Rand, array, randomNumber * sizeof(REAL), cudaMemcpyDeviceToHost) );

  FILE* out = Fopen(namefile, "w");

  for(int x=0;x<randomNumber;x++) { fprintf(out,"%f\n",h_Rand[x]); }

  fclose(out);

}


////////////////////////////////////////////////////////////////////////////////
// Write MT_RNG_COUNT vertical lanes of NPerRng random numbers to *d_Random.
// For coalesced global writes MT_RNG_COUNT should be a multiple of warp size.
// Initial states for each generator are the same, since the states are
// initialized from the global seed. In order to improve distribution properties
// on small NPerRng supply dedicated (local) seed to each twister.
// The local seeds, in their turn, can be extracted from global seed
// by means of any simple random number generator, like LCG.
////////////////////////////////////////////////////////////////////////////////
__global__ void RandomGPU( REAL *d_Random, int mode ){

  const int tid = blockDim.x * blockIdx.x + threadIdx.x;
  const int THREAD_N = blockDim.x * gridDim.x;

  int iState, iState1, iStateM, iOut;
  unsigned int mti, mti1, mtiM, x;
  unsigned int mt[MT_NN];

  REAL temp_rand;

  for(int iRng = tid; iRng < MT_RNG_COUNT; iRng += THREAD_N){
    //Load bit-vector Mersenne Twister parameters
    mt_struct_stripped config = ds_MT[iRng];

    //Initialize current state
    mt[0] = config.seed;
    for(iState = 1; iState < MT_NN; iState++)
      mt[iState] = (1812433253U * (mt[iState - 1] ^ (mt[iState - 1] >> 30)) + iState) & MT_WMASK;

    iState = 0;
    mti1 = mt[0];

    for(iOut = 0; iOut < NPerRng; iOut++){
      //iState1 = (iState +     1) % MT_NN
      //iStateM = (iState + MT_MM) % MT_NN
      iState1 = iState + 1;
      iStateM = iState + MT_MM;
      if(iState1 >= MT_NN) iState1 -= MT_NN;
      if(iStateM >= MT_NN) iStateM -= MT_NN;
      mti  = mti1;
      mti1 = mt[iState1];
      mtiM = mt[iStateM];

      x = (mti & MT_UMASK) | (mti1 & MT_LMASK);
      x = mtiM ^ (x >> 1) ^ ((x & 1) ? config.matrix_a : 0);
      mt[iState] = x;
      iState = iState1;

      //Tempering transformation
      x ^= (x >> MT_SHIFT0);
      x ^= (x << MT_SHIFTB) & config.mask_b;
      x ^= (x << MT_SHIFTC) & config.mask_c;
      x ^= (x >> MT_SHIFT1);

      temp_rand = ((REAL)x + 1.0f) / 4294967296.0f; //Convert to (0, 1] float

      if(mode==ZEROONE) { // write to global memory
        d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand;
      } else if(mode==ONEONE) {
        d_Random[iRng + iOut * MT_RNG_COUNT] = 2.0f*temp_rand-1.0f;
      }else {
        d_Random[iRng + iOut * MT_RNG_COUNT] = temp_rand - 0.5f;
      }

    }
  }
}


//-------------------------------------------------------------------------------------------------------
// HeatBath
//-------------------------------------------------------------------------------------------------------
#if defined(SINGLEHBKERNEL)
__global__ void HeatBath(REAL *a, REAL *b, REAL *c,
                         REAL *an, REAL *bn, REAL *cn,
                         REAL *Jpx, REAL *Jpy, REAL *Jpz,
                         REAL *Jmx, REAL *Jmy, REAL *Jmz,
                         REAL *atry, REAL *btry, REAL *ctry,
                         REAL *random, int even) {

  unsigned int iy, iz;

  int diff;
  REAL as, bs, cs;
  REAL anew, bnew, cnew;
  REAL hk, bhk, projection, pl, invnorm;

  REAL invexpbhk2;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x;


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    anew=atry[i];
    bnew=btry[i];
    pl=anew*anew+bnew*bnew;
    if(pl>d_const_zerotwentyfive) {
      continue;
    }
    cnew=ctry[i];
    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(((iz&1)^(iy&1))^even);

    spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) & d_maskL3;
    smz=(i-d_Lsdb2) & d_maskL3;

    as=an[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
      an[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
      an[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
      an[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
      an[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
      an[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];

    bs=bn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
      bn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
      bn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
      bn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
      bn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
      bn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];

    cs=cn[spx]*Jpx[(i&jindred[even])+(spx&jindblue[even])]+
      cn[spy]*Jpy[(i&jindred[even])+(spy&jindblue[even])]+
      cn[spz]*Jpz[(i&jindred[even])+(spz&jindblue[even])]+
      cn[smx]*Jmx[(i&jindred[even])+(smx&jindblue[even])]+
      cn[smy]*Jmy[(i&jindred[even])+(smy&jindblue[even])]+
      cn[smz]*Jmz[(i&jindred[even])+(smz&jindblue[even])];

    cs+=d_h;

    hk=SQRT(as*as+bs*bs+cs*cs);
    bhk=d_Beta*hk;

    invexpbhk2=EXP(-d_const_two*bhk);
    projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
    if(projection<-d_const_one) { projection=-d_const_one; }

    invnorm=SQRT((d_const_one-cnew*cnew)/pl);

    anew*=invnorm;
    bnew*=invnorm;

    factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);

    anew = anew-factor*as;
    bnew = bnew-factor*bs;
    cnew = cnew-factor*cs;

    factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
    projection/=hk;
    a[i]=as*projection+anew*factor;
    b[i]=bs*projection+bnew*factor;
    c[i]=cs*projection+cnew*factor;

  }

}

#else

//-------------------------------------------------------------------------------------------------------
// HeatBathRed
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathRed(REAL *a, REAL *b, REAL *c,
                            REAL *an, REAL *bn, REAL *cn,
                            REAL *Jpx, REAL *Jpy, REAL *Jpz,
                            REAL *Jmx, REAL *Jmy, REAL *Jmz,
                            REAL *atry, REAL *btry, REAL *ctry,
                            REAL *random) {

  unsigned int iy, iz;

  int diff;
  REAL as, bs, cs;
  REAL anew, bnew, cnew;
  REAL hk, bhk, projection, pl, invnorm;

  REAL invexpbhk2;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x;


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    anew=atry[i];
    bnew=btry[i];
    pl=anew*anew+bnew*bnew;
    if(pl>d_const_zerotwentyfive) {
      continue;
    }
    cnew=ctry[i];
    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    diff=(iz&1)^(iy&1);

    spx=(i >> (-diff)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (diff-1));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (-diff)) + (i>>(diff-1));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) & d_maskL3;
    smz=(i-d_Lsdb2) & d_maskL3;

    as=an[spx]*Jpx[i]+
      an[spy]*Jpy[i]+
      an[spz]*Jpz[i]+
      an[smx]*Jmx[i]+
      an[smy]*Jmy[i]+
      an[smz]*Jmz[i];

    bs=bn[spx]*Jpx[i]+
      bn[spy]*Jpy[i]+
      bn[spz]*Jpz[i]+
      bn[smx]*Jmx[i]+
      bn[smy]*Jmy[i]+
      bn[smz]*Jmz[i];

    cs=cn[spx]*Jpx[i]+
      cn[spy]*Jpy[i]+
      cn[spz]*Jpz[i]+
      cn[smx]*Jmx[i]+
      cn[smy]*Jmy[i]+
      cn[smz]*Jmz[i];

    cs+=d_h;

    hk=SQRT(as*as+bs*bs+cs*cs);
    bhk=d_Beta*hk;

    invexpbhk2=EXP(-d_const_two*bhk);
    projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
    if(projection<-d_const_one) { projection=-d_const_one; }

    invnorm=SQRT((d_const_one-cnew*cnew)/pl);

    anew*=invnorm;
    bnew*=invnorm;

    factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);

    anew = anew-factor*as;
    bnew = bnew-factor*bs;
    cnew = cnew-factor*cs;

    factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
    projection/=hk;
    a[i]=as*projection+anew*factor;
    b[i]=bs*projection+bnew*factor;
    c[i]=cs*projection+cnew*factor;

  }

}


//-------------------------------------------------------------------------------------------------------
// HeatBathBlue
//-------------------------------------------------------------------------------------------------------
__global__ void HeatBathBlue(REAL *a, REAL *b, REAL *c,
                             REAL *an, REAL *bn, REAL *cn,
                             REAL *Jpx, REAL *Jpy, REAL *Jpz,
                             REAL *Jmx, REAL *Jmy, REAL *Jmz,
                             REAL *atry, REAL *btry, REAL *ctry,
                             REAL *random) {

  unsigned int iy, iz;

  int diff;
  REAL as, bs, cs;
  REAL anew, bnew, cnew;
  REAL hk, bhk, projection, pl, invnorm;

  REAL invexpbhk2;

  int spx, smx, spy, smy, spz, smz;

  REAL factor;

  int i;

  const unsigned int tid = threadIdx.x;

  const unsigned int num_threads = gridDim.x*blockDim.x;
  int istart;

  istart=blockDim.x*blockIdx.x;


  for(i=istart+tid; i<d_nspin; i+=num_threads) {

    anew=atry[i];
    bnew=btry[i];
    pl=anew*anew+bnew*bnew;
    if(pl>d_const_zerotwentyfive) {
      continue;
    }
    cnew=ctry[i];
    iz=i>>d_sdbLsm1;
    iy=(i - (iz<<d_sdbLsm1))>>d_sdbLm1;

    diff=(iz&1)^(iy&1);

    spx=(i >> (diff-1)) + (((i-(i&d_maskL)) + ((i+1)&d_maskL)) >> (-diff));
    smx=(((i-(i&d_maskL)) + ((i-1)&d_maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&d_maskL2))+((i+d_Ldb2)&d_maskL2);
    smy=(i-(i&d_maskL2))+((i-d_Ldb2)&d_maskL2);
    spz=(i+d_Lsdb2) & d_maskL3;
    smz=(i-d_Lsdb2) & d_maskL3;

    as=an[spx]*Jpx[spx]+
      an[spy]*Jpy[spy]+
      an[spz]*Jpz[spz]+
      an[smx]*Jmx[smx]+
      an[smy]*Jmy[smy]+
      an[smz]*Jmz[smz];

    bs=bn[spx]*Jpx[spx]+
      bn[spy]*Jpy[spy]+
      bn[spz]*Jpz[spz]+
      bn[smx]*Jmx[smx]+
      bn[smy]*Jmy[smy]+
      bn[smz]*Jmz[smz];


    cs=cn[spx]*Jpx[spx]+
      cn[spy]*Jpy[spy]+
      cn[spz]*Jpz[spz]+
      cn[smx]*Jmx[smx]+
      cn[smy]*Jmy[smy]+
      cn[smz]*Jmz[smz];

    cs+=d_h;

    hk=SQRT(as*as+bs*bs+cs*cs);
    bhk=d_Beta*hk;

    invexpbhk2=EXP(-d_const_two*bhk);
    projection=d_const_one+(LOG(invexpbhk2+(d_const_one-invexpbhk2)*random[i]))/bhk;
    if(projection<-d_const_one) { projection=-d_const_one; }

    invnorm=SQRT((d_const_one-cnew*cnew)/pl);

    anew*=invnorm;
    bnew*=invnorm;

    factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);

    anew = anew-factor*as;
    bnew = bnew-factor*bs;
    cnew = cnew-factor*cs;

    factor=SQRT((d_const_one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew));
    projection/=hk;
    a[i]=as*projection+anew*factor;
    b[i]=bs*projection+bnew*factor;
    c[i]=cs*projection+cnew*factor;

  }

}


#endif
